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ABSTRACT 

This paper provides a simple, efficient, and robust numerical technique 
for solving two-dimensional incompressible steady viscous flows at moderate- 
to-high Reynolds numbers. The proposed approach employs an incremental 
multigrid method and an extrapolation procedure based on minimum residual 
concepts to accelerate the convergence rate of a robust block-line-Gauss- 
Seidel solver for the vorticity-stream function Navier-Stokes equations. 

Results are presented for the driven cavity flow problem using uniform 
and nonuniform grids and for the flow past a backward facing step in a 
channel. For this second problem, mesh refinement and Richardson 
extrapolation are used to obtain useful benchmark solutions in the full range 
of Reynolds numbers at which steady laminar flow is established. 
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INTRODUCTION 


This paper is concerned with the simulation of two-dimensional 
incompressible steady laminar separated flows at moderate-to-high Reynolds 
numbers (Re), using a simple, efficient, and robust numerical technique. 
Among the many numerical methods developed for the incompressible Navier- 
Stokes equations, those recently employed to solve high Re steady separated 
flows are very complex and sophisticated. For example, (i) Ghia et al. [1] 
use the cumbersome coupled strongly implicit method as a robust smoother for 
the already involved full-approximation-storage, full-raultigrid method of 
Brandt [2], and (ii) Schreiber and Keller [3] solve a fourth order nonlinear 
problem for the stream function by a sequence of Newton and chord iterations, 
and use a costly L-U factorization with partial pivoting to solve the large 
sparse linear systems associated with the Newton iteration. In both 
techniques, the solution at a lower value of Re is to be used effectively to 
generate a sufficiently good initial condition. Therefore, it appears 
worthwhile to provide a numerical technique for solving high-Re separated 
flows, which is possibly as powerful and efficient as the best methods 
available to date but much simpler to implement and to use. 

In the last few years, the second author has developed approximate 
factorization [4] and line relaxation [5] methods for solving the steady-state 
vorticity-stream function Navier-Stokes equations. These methods employ a 
two-level implicit Euler time stepping and the delta form [6] to discretize 
and linearize in time the unsteady governing equations and make effective use 
of a deferred correction strategy for the finite-difference spatial 
discretization; namely, second-order-accurate central differences are used for 
all spatial derivatives except the advection terms in the left hand side (LHS) 
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iraplicit operator, which are discretized using first-order-accurate upwind 
differences. In this way, an artifical viscosity is introduced which is 
proportional to a time derivative and thus vanishes as the sought steady-state 
solution is reached (see Appendix A). Also, the large 2x2 block-pentadiagonal 
matrix associated with the LHS implicit operator is diagonally dominant, so 
that the Alternating Direction Implicit (ADI) [4] or line Gauss-Seidel (LGS) 
[5] solution procedures enjoy the robustness and stability of upwind schemes 
and the accuracy of central-difference schemes. Both methods are very simple 
and have been reasonably successful in computing steady flows at moderate 
Reynolds numbers. However, their convergence rate invariably deteriorates 
when the computational mesh is refined and/or the Reynolds number is 
increased. In an attempt to overcome such a limitation, an incremental 
multigrid approach has been recently proposed [7], which is particularly 
suitable for this type of numerical methods, extremely simple, and does not 
require any additional storage with respect to the basic numerical scheme used 
as a smoother, nor any sophisticated strategy for cycling among the various 
grids. Therefore, it could be a viable alternative to more complicated 
multilevel methods. However, its validity has only been demonstrated for a 
model problem and is restricted to the case of uniform grids. It seems 
therefore necessary and appropriate to assess its merits and deficiencies 
versus more difficult problems and to further improve its performance, without 
affecting its major merit, namely, its simplicity. 

These goals are achieved in this paper, which: (i) provides an improved 
version of the incremental multigrid method of [7], capable of handling meshes 
with reasonably high stretching; (ii) supplements such a method with an 
extrapolation technique based on minimum residual concepts [8] to further 
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enhance its efficiency; (iii) employs the resulting procedure to provide a 
benchmark solution for flow past a backward facing step in a channel in the 
full range of Re at which steady laminar flow is established. 


NUMERICAL METHOD 

The nondimens ional vorticity-stream function Navier-Stokes equations are 
given in the standard Cartesian coordinate system, for simplicity, as 


<l). + ill (i) -ill 0) - ■=— (u +u ) = 0 
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In Eqs. (1-2), Re is the Reynolds number, id and i(i are the vorticity and 

the stream function, t is the time, x and y are the horizontal and 

vertical Cartesian coordinates, and subscripts indicate partial derivatives. 

Equations (1-2) are discretized in time by means of a two level implicit Euler 

time stepping and linearized using the delta approach [6], by neglecting terms 
2 

of order A , to give: 
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where At is the time step, the superscript n indicates the known 

solution at the time level t n and Aw, Aip are the unknowns to be 

computed. Equations (3-4) are discretized in space using second-order- 

accurate central differences throughout, except for the advection 

derivatives Aw , A^ , Aw and Aib in the LHS of Eq. (3) which are 

discretized using first-order-accurate upwind differences according to the 

signs of ip n , w n , i|i n , w n , and solved approximately by a block-ADI [4] or 
y x x y 

block-LGS method [5]. Only block-tridiagonal systems need to be solved along 
each row and column of the computational domain, and the double boundary 
condition for \jj can be easily imposed to provide the value of the 

vorticity at the wall directly (see Appendix B). Two points are of 
interest: (i) a relaxation-like time derivative needs to be added to the 

stream function equation if an ADI solution procedure is employed [9]; (ii) 
the advection terms in the right hand side of Eq. (4) are replaced by the 


corresponding conservative form 


- (♦ w) + (4> a>) (5) 

y x x y 

which has been shown to provide more accurate results (see, e.g., [5]). This 
amounts again to employing a deferred correction approach, which is made 
particularly elegant and simple to implement by the use of the delta form 
[4]. Notice, in fact, that a standard central difference discretization of 
Eq. (5) requires values of \p from the NW (North-West), SE (South-East), NE 
and SW gridpoints in the computational stencil and, if used in the implicit 
left hand side operator, would increase the number of nonzero diagonals in the 
resulting matrix and reduce its diagonal dominance. After every ADI or LGS 
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sweep, the solution is advanced as 

(ij/ , w) n *■ (i|», w) n + (Ai|> , Aw) 


( 6 ) 


and the process is repeated until a satisfactory convergence criterion is met. 

In order to describe the multigrid procedure employed in this study, Eqs. 
(3-4) are rewritten in a more general form, by dropping the superscript n 
and introducing superscripts H and h to indicate the current and the 
finest grids used in the computations (H = h, 2h, 4h, and 8h) 
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where indicates the standard 9-point collection operator, applied as 

many times as needed to go from the finest mesh h to the current mesh H. 

Starting from an arbitrary initial condition, Eqs. (7-8) are solved on the 

finest grid h — where they coincide with Eqs. (3-4) — by means of a two 

ti h 

sweep alternating direction block-LGS iteration, to provide Aw , Atp ; the 
Vi Vi 

solution w , f is updated and Eqs. (7-8) are solved on successively 
coarser grids (H = 2h, 4h, and 8h); the entire process is repeated until the 
finest-grid residual is reduced to a suitably small value. In more detail, at 
every grid level H, the following steps are required by the proposed 
multigrid strategy: a) the coefficients in the LHS of Eqs. (7-8) are 
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Vi ti 

evaluated at the H-mesh gridpoints using the finest-grid solution (to , ) 

locally, whereas the RHS steady state residuals are evaluated on the finest 

grid h and collected up to the current grid H; b) Eqs. (7-8) are then 

solved approximately, using a single sweep of the aforementioned block-LGS 

H H 

smoother and homogeneous Dirichlet boundary conditions, to provide Ato , A^i ; 
tl tl 

c) Aw , Ai|) are evaluated as 

(Ao) h , Ai|> h ) = ijj (Aoj H , A<|> H ) (9) 

where ijj is the standard bilinear interpolation operator from the current 
grid H to the finest grid h; d) the finest-grid solution is updated as 

(w h , ip* 1 ) «- (w h , *|> h ) + (Am h , A^ h ); (10) 

e) the vorticity at the boundaries is finally corrected so as to satisfy the 

no-slip boundary condition on the finest mesh (see Appendix B). All of the 

aforementioned steps are performed twice, with the block-LGS solution method 
marching from left to right and from top to bottom of the computational 
domain, respectively. A multigrid cycle is shown schematically in Figure 1, 

where it is seen to differ from both the more usual V and saw-tooth 
cycles. It is noteworthy that the proposed methodology is very simple, since 
it does not require any logical choices to be made and employs a single free 
parameter, namely, the time step At. Furthermore, it does not need any 
additional storage with respect to the basic smoother, insofar as only the 
finest-grid solution is computed and a single array is used for the deltas at 
all grid levels. However, its work per iteration is slightly greater than 


that required by most current multigrid methods, due to the additional 
interpolations and collections needed to visit and update the finest-grid 
solution after every coarse-grid calculation and, due to its extreme 
simplicity, it is likely to be less efficient than more sophisticated 
multigrid methods. 

The present approach, as described above, can be applied without any 
modifications to the vorticity-stream function equations written in a general 
curvilinear coordinate system £, q. The scale factors and the Jacobian of 
the transformation (x, y) ♦ (g , q) are evaluated once and for all on the 
finest mesh and treated as the other variable coefficients in the linearized 
discrete equations arising at every grid level. However, numerical 
experiments performed for the case of the driven cavity flow problem have 
shown that the efficiency of the method rapidly deteriorates as the 
computational grid in the physical plane becomes increasingly nonuniform. 
Therefore, following the lead of several other workers (see, e.g., [10]), the 
9-point collection operator for the residual has been modified so as to use 
weighed areas in physical space, and the bilinear interpolation operator has 
been modified so as to use distances among gridpoints also in physical 
space. More precisely, in order to collect a quantity f from the finest 

mesh h to the mesh H = 2h at point P, the standard 9-point collection 
operator is given as: 


C? h f {4f 4 . + 2f . . . + 2f, . .. + 2f . . + 2f . . . + 

h 16 1 i,j l-l, j i,j+l i+l, J i»J-l 


( 11 ) 
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whereas the modified collection operator is 
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where Aj, A 2 , A^ and A^ are the areas of the four cells surrounding the 
gridpoint P (see Figure 2 which shows the 9-point computational stencil in 
the physical (x, y) and computational (£ , p) planes). On the other hand, 
in order to interpolate a quantity f at point Q using the ^i j an< * 
^i+2 j values available on the coarser mesh, the standard bilinear 
interpolation operator is given as 


f + f 

h f _ 1 i+2 , j r i,j 

2h 1 2 


(13) 


whereas the modified interpolation operator is 



f i+2,j (x i+l ~ x i> + f i, j (x i+2 ~ x i+l> 
x i+2 ~ X i 


(14) 


Finally, in order to further enhance the convergence rate of the method, 
the following extrapolation technique based on minimum residual concepts [8] 
is used, after every k multigrid cycles, to obtain a new initial condition 
for the finest-mesh solution. Let fn-2^ f n be the solution vectors 
(the vectors of all co and ^ gridpoint values) at the end of the last 
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three cycles and R n 2 , R n 1 > R n the corresponding residuals. A new 

* 

initial solution f is obtained as 


f* - £ n -‘ + Cl (f" - f"’ 1 ) + - f"' 2 ) 


with £ ^ and ^ evaluated as follows. The residual R is assumed to 


depend linearly on and ? 2> as 


R* - R n 1 + ^(r” - R n L ) + C 2 (R n-1 - R n 2 ) 


and the dot product 


•k k 

R • R is minimized with respect to and ? 2 


to give: 


a b\ /c 


b c/ \ c 


where 


a = (R n - R n_1 ) • (R n - R n_1 ) 


b = (R n 1 - R n ~ 2 ) • (R n - R n l ) 


c = (R n_1 - R n 2 ) • (R n 1 - R n 2 ) 


d - R n 1 • (R n - R n_1 ) 


e = R n_1 • (R n_1 - R n “ 2 ). 


( 22 ) 
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It is noteworthy that such a procedure, which can be implemented using an 
arbitrary number of extrapolation parameters [8], is extremely simple 
and employs a negligible amount of CPU time with respect to the basic 
solver* On the other hand, it requires additional memory, insofar as both the 
solution and the residual vectors are needed at previous iteration levels, and 
introduces an additional parameter in the proposed numerical method, namely, 
the interval of application of the extrapolation procedure, k. However, 
memory is not a problem, especially for the present case of two-dimensional 
flows and the convergence rate of the method has been found here to be rather 
insensitive to the value of k (see also [8]). A final remark is needed* In 
the present study, both the two-parameter extrapolation described above and 
the simpler one based on a single parameter 5^ have been employed* The 
two-parameter technique has consistently provided better results, but the 
efficiency gain achieved with respect to the simpler one-parameter approach 
has been rather limited, so that no attempt at using three or more parameters 
was made. 


RESULTS 

The numerical technique, as described in the previous section, has been 
applied to solve two viscous flow problems for several values of Re. The 
computations were always started from rest and used a nonoptimized time step, 
usually equal to one. 
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Flow in a driven cavity 

The classical driven cavity flow [11] was considered at first in order to 
assess the performance: (i) of the basic multigrid method for increasing 
values of Re, without and with the extrapolation procedure; (ii) and of the 
modified method for increasingly nonuniform grids. As far as the first point 
is concerned, calculations were performed for Re = 1000 using a 129x129 
uniform grid and from one to four grid levels, without and with the 
extrapolation technique applied every 20 iterations. The convergence 
histories are given in Figures 3 and 4, where the logarithm of the (Lq norm 
of the) vorticity residual is plotted versus the work units, one work unit 
being the CPU time required to complete a two-sweep iteration on the finest 
mesh. In all cases, the residual has been dropped to machine zero on a Gould 
PN9005 computer using single precision arithmetic. It clearly appears that 
the multigrid method provides a considerable improvement over the basic 
smoother and that the extrapolation technique further enhances its 
efficiency. In order to assess the influence of the interval of application 
of the extrapolation procedure, k, on the convergence rate of the method, 
results have been obtained for various values of k and are given in Table 1, 
as the work units necessary for the vorticity residual to reach 10 The 
value of k is seen to have a minor influence on the convergence rate of the 
method (see also [8]) and can thus be chosen somewhat arbitrarily. 
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k 

10 

15 

20 

25 

work units 

149 

167 

145 

161 


Table 1. Influence of k on the convergence rate of the method 

The more difficult Re = 3200 flow case was then considered in order to 
further test the robustness of the method. Figure 5 provides the convergence 
histories of the basic smoother and of the four-grid multigrid method without 
and with the extrapolation procedure applied every 20 iterations. The basic 
smoother, although stable, is extremely slow to converge and also the 
raultigrid method experiences rather severe difficulties before being able to 
reduce the residual effectively. Also, due to the lack of smoothness in the 
convergence history of the scheme, the extrapolation procedure is found to 
actually delay convergence. Incidentally, for Re = 10,000, convergence 
requires more than 10,000 work units, the extrapolation procedure again being 
beneficial. In conclusion, the present multigrid method is extremely robust 
but becomes inefficient for very high values of Re. 

In order to address the second point of interest, namely, the performance 
of the improved method for the case of nonuniform grids, the same driven 
cavity problem was considered, again for Re « 1000 and 3200, by mapping the 
physical plane into a uniform-grid computational domain using the following 
analytical transformation [4], [5]: 
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(*) = 0.5 + 0.5 tanh [C(^ I })]/tanh(C). (23) 

For C << 1, the x and y lines are practically equally spaced, whereas, 

as C increases, more and more gridlines are concentrated near the boundaries 

of the unit-square physical domain. The governing equations in terms of the £ 

and n variables are given in [5] and the scale factors and the Jacobian of 

the transformation (23) are computed numerically using second-order-accurate 

central differences everywhere except at the boundaries, where three-point 

one-sided differences are used [5]* In the present calculations a 65x65 

uniform grid in the £ , n computational plane was used, for several values 

of C, and a reduced value of the time step, At » 0.2, was always employed, 

as already in [5]. The improved four-grid raultigrid method converged without 

any difficulty for C as high as 1.4, for which the maximum- to-minimura 

Ax (Ay) ratio is equal to 4.45. Also, the extrapolation procedure improved 

the efficiency of the method for both Re = 1000 and Re = 3200, convergence 

to machine zero requiring about 400 and 1000 work units, respectively. For 

completeness, the numerical results are given in Table 2 as the maximum values 

of the stream function and the values of the vorticity at the center 

M 

of the moving plate (w p ). The corresponding results obtained using 

uniform grids of 97x97 and 129x129 gridpoints are also given for 

comparison. The 65x65 nonunif orm-grid results are as accurate as the 

129x129 uniform-grid ones, so that, for the present problem, the nonunif orra- 
grid method results to be more effective overall. However, the present 

approach is considered to be inadequate to compute external flows requiring 
highly stretched grids. 
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W C 


97x97 

uniform 

.1174 

14.95 

Re = 1000 

129x129 

uniform 

.1180 

14.88 


65x65 

nonuniform 

.1181 

14.88 



97x97 

uniform 

.1166 

26.98 

Re = 3200 

129x129 

uniform 

.1187 

26.16 


65x65 

nonuniform 

.1193 

25.96 


Table 2* Driven cavity flow results 


Flow past a backward facing step 

The flow past a backward facing step in a channel, see Figure 6, is a 
very interesting problem which has been chosen by the organizers of a GAMM 
workshop as the test case for comparing a great number of codes for solving 
the incompressible Navier-Stokes equations. From the results presented at the 
workshop [12] , it clearly appears that for Re >_ 500 most methods face 
convergence difficulties and/or need some kind of upwinding to handle the flow 
regions where convection dominates diffusion. Physically, as clearly shown by 
the very careful experiments of Armaly et al. [13], the structure of the flow 
becomes more and more complicated as Re increases: the flow, which always 
separates over the step, reattaches downstream at a distance which increases 
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with Re and, for sufficiently high values of Re, a secondary separation 
region develops on the wall opposite to the step. This problem was thus 
chosen as a severe test for the present approach using central differences for 
the RHS steady state residual. The computational domain is limited to the 
interior of the channel immediately at the right of the step and a fully 
developed (Couette flow) parabolic velocity profile is used as a boundary 
condition in the upper half of the left boundary (h/H = 0.5, see Figure 6) 
[13, 14]. The proposed approach, without and with the extrapolation technique 
has been employed using uniform grids with 49x49, 65x49, 81x49 and 97x49 
gridpoints for the cases Re = 200, 400, 600 and 800, respectively, the 
downstream boundary condition being set at a distance from the step equal to 
7.5, 10, 12.5 and 15. The nondimensional height of the channel H is equal 
to 1 and the maximum value of the nondimensional longitudinal velocity 
component at inlet is equal to 1.5 [14], At the outlet of the channel, 

second-order-accurate three-point homogeneous Neumann boundary conditions are 
used for both ^ and u), to minimize the upstream influence due to imposing 
an asymptotic condition at a finite downstream distance. At the inlet and at 
all of the walls, standard no-slip conditions are prescribed, as shown in 
Appendix B. In particular, at the inlet, ’I'yy. which is discontinuous at the 
corner C, is evaluated analytically, with the gridpoint C obviously being 
considered part of the inlet-flow domain ^yy = cases > no 

convergence difficulty was encountered, again starting all computations from 
rest and always using At = 1. The convergence histories for the method, 
using from 1 to 4 grid levels, without and with the extrapolation applied 
every 20 iterations, are given in Figures 7 and 8 for Re = 200, and in 
Figures 9 and 10 for Re = 800. For the simpler Re = 200 flow case, using a 
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rather coarse 49x49 finest mesh, the multigrid approach reaches its peak 
efficiency when using 3 grid levels, without the extrapolation, and 2 grid 
levels, with the extrapolation. For the more difficult Re = 800 flow case, 
using a 97x49 finest-grid, the efficiency of the multigrid method always 
improves with increasing number of grid levels. In all cases, the 
extrapolation significantly improves the performance of the approach. 

An efficient and second-order-accurate method being available, solutions 
were obtained for all four cases doubling the number of mesh intervals in both 
directions, so as to provide a benchmark solution for this very interesting 
problem. Figures 11 and 12 show the lower and upper walls vorticity 
distributions obtained using 97x97, 129x97, 161x97 and 193x97 gridpoints 
for Re = 200, 400, 600 and 800, respectively. On the same figures, the 
results obtained using the coarser grids are also given as symbols. It 
appears that for Re = 200 and 400 grid convergence has been achieved, 
whereas for Re = 600 and 800 further mesh refinement is probably 
warranted. However, the two different grid results in Figures 11 and 12 are 
reasonably close, so that Richardson extrapolation can be used with confidence 
to obtain a benchmark solution: Table 3 provides the values of the locations 
of the reattachment point for the primary separation bubble (X1R) and of the 
separation and reattachment of the secondary separation bubble (X2S, X2R) , 
divided by the height of the step h [14], obtained using linear interpolation 
between the two gridpoints at which the wall vorticity changes sign and 
Richardson extrapolation to zero step size. Incidentally, the numerical 
results used for the extrapolation are converged to machine zero, using double 
precision arithmetic. 
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Re 

200 

400 

600 

800 

XIR/h 

5.34 

8.63 

10.72 

12.16 

X2S/h 

- 

7.96 

8.71 

9.67 

X2R/h 

- 

10.44 

16.23 

20.96 


Table 3. Benchmark Results 

It needs to be remarked that, for all values of Re, the far downstream 
values of the vorticity on the lower and upper walls should be 3 and -3, 
respectively. From the results of Figures 11 and 12, one may thus believe 
that the outflow boundary conditions have not been imposed far enough 
downstream, especially for the higher values of Re. Therefore, the coarser 
grid computations were repeated for the cases Re = 200 and Re ■ 800, moving 
the outflow boundary-condition locations to x = 15 and x = 25, 
respectively, and increasing the number of longitudinal gridpoints to maintain 
the same value of Ax. The results for the lower and upper walls vorticity 
are given in Figures 13 and 14 for both sets of calculations. The vorticity 
is seen to tend to its asymptotic value correctly and the results obtained 
using the two different locations for the outflow boundary conditions are in 
perfect agreement. The usefulness of using outflow conditions of Neumann type 
is thus clearly demonstrated so as the validity of the results in Table 3 as a 


benchmark solution. 
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CONCLUSIONS 

A simple, robust, and efficient method has been developed for solving 
two-dimensional steady viscous flows. An incremental alternating direction, 
block-line-Gauss-Seidel relaxation method using first-order-accurate upwind 
differences in the left hand side implicit operator and second-order-accurate 
central differences in the right-hand-side steady-state residual is used as 
smoother within a very simple multigrid algorithm, supplemented by an 
extrapolation procedure based on minimum residual concepts. The proposed 
technique has been tested versus the classical driven cavity flow, for values 
of the Reynolds number (Re) as high as 10,000, and used to provide useful 
benchmark solutions for flow past a backward facing step in a channel, for 
values of Re covering the full range at which steady laminar flow exists. The 
convergence rate of the method, which always starts from an arbitrary initial 
condition and marches towards steady state using a simple multigrid cycle 
without any optimization, logical choices or adjustable parameters, is very 
satisfactory for moderate-to-high values of Re. However, a more 
sophisticated approach is required for very high values of Re and/or highly 
nonuniform grids. 
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APPENDIX A 


Consider the linear advection diffusion equation 


u + cu - eu = 0 
t x xx 


(Al) 


where c is a function of x and can be either positive or negative. The 
discrete form of Eq. (Al) using the delta approach and a deferred correction 
strategy as done in this paper for the vorticity-stream function equations is 


Au^ ^ c i + |cj A Uj[ - Au i _ 1 ^ c t - IcJ 
At ’ 2 Ax 2 


Au 


i+1 


- Au . 


Ax 


(A2) 


- e 


au i+i ■ “ u i + au i-i 


u 


= “C, 


i+1 


U i-1_ + e U i+1 2u i+l u i-l 


2Ax 


Ax 


The two incremental advection terms in Eq* (A2) can be written as 


Au i+1 " Au i-1 
C i IKx 


c^ Ax Au.^ 
2 


2Au^ + Au^ + j 


Ax 


(A3) 


so that Eq. (A2) is easily seen to be an implicit central-in-space finite 
difference discretization of Eq. (Al), plus an artificial viscosity term which 
is the backward-in- time central-in-space finite difference approximation of 


- | c i | Ax At 

u 

xxt 


(A4) 


2 
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and thus vanishes identically at steady state. Similarly, it is seen that the 
four advective terms in the LHS of Eq. (3) are equivalent to the corresponding 
central difference approximations plus artificial viscosity terms which vanish 
at steady state. It is to be pointed out that the discretization used in Eq. 
(A2) is the delta form of the one proposed by Khosla and Rubin [16] and is 
easily seen to provide a diagonally dominant matrix for the LHS implicit 
operator. 


-24- 


APPENDIX B 


Let us consider the gridpoints adjacent to the boundary line BB, together 
with a mirror image point, 0, outside the computational domain [11] as shown 
in Figure 15. At gridpoint 1, the double specification for the stream 

function is given as: 

'J'j = a (Bl) 

('[' x ) j = b. (B2) 

Equation (B2) is discretized using a third-order-accurate four point 

difference 


(- 4)3 + 6 tj>2 “ 3 ^1 - 2 ^q) 

g__ 


= b. 


(B3) 


In order to eliminate the additional unknown (ji^, the steady state stream 
funtion equation is also used at the boundary gridpoint 1 [11] 


\b + + u) = 0 

T W T 


XX 


yy 


(B4) 


which is discretized as: 


'J'o " 2^1 + V 2 


Ax 


+ t|> + oj , = 0 

yy 1 


(B5) 


where the iji term is left unchanged for convenience. By combining Eqs. 

(B3) and (B5), the following equation for the vorticity at the boundary, oi^, 
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is obtained: 


a) , = — ib + 

i yy 


6 b Ax + 7 ip j — 8 2 'I' 3 


2 Ax 


(B6) 


Equations (Bl) and (B6) are written in delta form and used, together with the 
internal-gridpoints discrete equations and the corresponding conditions for 
the RHS boundary, to provide a 2 x 2 block-tridiagonal system which is solved 
very efficiently by block-tridiagonal elimination. Notice that in Eq. (B6) 

^|i is either zero, if line BB is a solid boundary, or is known, if line BB 
is a flow-inlet boundary. Also, from Eq. (B6), it clearly appears that a 
third order accurate discretization of Eq. (B 2 ) is needed to obtain a second 
order accurate o>^ (see also [ 1 ]). Finally, in the present multigrid 
method, Eq. (B6) and the corresponding ones are also used to correct the 
finest-grid solution at the boundaries, after every coarse-to-f ine-grid 


interpolation. 
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Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 


FIGURE CAPTIONS 

1. Schematic of the multigrid cycle. 

2. Finite difference stencil in physical (x,y) and 
computational (5 ,n ) planes. 

3. Convergence histories of the multigrid method using 1, 2, 3 and 
4 grid levels for Re = 1000. 

4. Convergence histories of the multigrid method with extrapolation 

using 1, 2, 3 and 4 grid levels for Re = 1000. 

5. Convergence histories of the basic solver and of the four-grid 

multigrid without and with extrapolation (dotted line) for Re = 
3200. 

6. Flow past a backward facing step in a channel: geometry and 

boundary conditions. 

7. Convergence histories of the multigrid method using 1, 2, 3 and 
4 grid levels for Re = 200. 

8. Convergence histories of the multigrid method with extrapolation 
using 1, 2, 3 and 4 grid levels for Re = 200. 

9. Convergence histories of the multigrid method using 1, 2, 3 and 
4 grid levels for Re = 800. 

10. Convergence histories of the multigrid method with extrapolation 
using 1, 2, 3 and 4 grid levels for Re = 800. 

11. Effect of grid refinement on the lower wall vorticity for 
various values of Re. 
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Flgure 12. Effect of grid refinement on the upper wall vorticity for 
various values of Re. 

Figure 13. Effect of downstream boundary condition location on the lower 
wall vorticity for two values of Re. 

Figure 14. Effect of downstream boundary condition location on the upper 
wall vorticity for two values of Re. 

Figure 15. Computational gridpoints around a boundary line. 
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